##############################
####
#### Authors: Stefan Müller, and Tom Louwerse
#### Title: "The Electoral Cycle Effect in Parliamentary Democracies."
#### Journal: Political Science Research and Methods
#### 
#### Replication Material
#### 
##############################


### File description: This file reproduces all regressions, plots, 
### and tables based on the nlme package.
### Note that the NLME models can take a long time to run. 


# 1) Load libraries and prepare data -----
library(tidyverse)
library(lme4)
library(effects)
library(texreg)
library(car)
library(pastecs)
library(stargazer)
library(ggthemes)
library(nlme)

# Load functions for texreg tables and ggplot2 theme
source("helper.R")

# Set theme for all plots
theme_set(theme_base_edited())

### This file runs all regressions based on the NLME package. 
### Note that it will take long to run these models.

# Load data
dta <- readRDS("data.rds")

# Apply filters
dta_all <- dta %>%
    dplyr::filter(total_length >= 182.2) %>%
    dplyr::filter(electoral_cycle_cabinet > 0, electoral_cycle_max < 1.2) %>%
    dplyr::filter(election == "Legislative") %>%
    dplyr::filter(system != "Presidential")

# Only select polls by government parties
dta_gov <- dta_all %>%
    dplyr::filter(government == 1 & electoral_cycle_cabinet > 0)


# Removing missing values required to plot resid vs fitted
dta_gov_base <- dta_gov %>% 
    filter(!is.na(poll_change) & !is.na(electoral_cycle_cabinet)
           & !is.na(gdp_change_lag) & !is.na(vote_share_last) & !is.na(elecyear))


#### NLME models -----

reg_base_nlme <- lme(fixed = poll_change ~ poly(electoral_cycle_cabinet, 3) +
                         gdp_change_lag +
                         vote_share_last +
                         I(elecyear - 1986),
                     random = ~ 1| country / election_id / partyid,
                     correlation = corCAR1(value = 0.2,
                                           form =  ~ electoral_cycle_cabinet | country / election_id / partyid),
                     data = dta_gov_base)

reg_base_nlme2 <- lme(fixed = poll_change ~ poly(electoral_cycle_cabinet, 2) +
                          gdp_change_lag +
                          vote_share_last +
                          I(elecyear - 1986),
                      random = ~ 1| country / election_id / partyid,
                      correlation = corCAR1(value = 0.2,
                                            form =  ~ electoral_cycle_cabinet | country / election_id / partyid),
                      data = dta_gov_base)

anova(reg_base_nlme, reg_base_nlme2)

# Prime minister dissolution power (NLME)
reg_pm_diss_nlme <- update(reg_base_nlme2, . ~ . + poly(electoral_cycle_cabinet, 2) * diss_pm,
                           control = lmeControl(msVerbose = TRUE, sing.tol = 1e-20))

## Single-party cabinet (NLME)
reg_single_party_gov_nlme <- update(reg_base_nlme2, . ~ . + 
                                        poly(electoral_cycle_cabinet, 2) * single_party_cabinet,
                                    control = lmeControl(msVerbose = TRUE, sing.tol=1e-20))


## Model decades (NLME)

# Select countries where we have polling data since the 1960s
country_selection <- c("Australia", "Canada", 
                       "Denmark", "Germany", 
                       "Netherlands", "Norway", 
                       "United Kingdom")

dta_decades <- subset(dta_gov, elecyear > 1959 & country_code %in% country_selection)

dta_decades <- dta_decades %>%
    dplyr::arrange(country, elecdecade) %>%
    filter(elecyear > "1959")

dta_decades <- dta_decades %>%
    filter(!is.na(poll_change) & !is.na(electoral_cycle_cabinet)
           & !is.na(gdp_change_lag) & !is.na(vote_share_last) & !is.na(elecyear))


model_decades_nlme <- lme(fixed = poll_change ~ 
                              poly(electoral_cycle_cabinet, 3) * elecdecade +
                              gdp_change_lag +
                              vote_share_last,
                          random = ~ 1| country / election_id / partyid,
                          correlation = corCAR1(form =  ~ electoral_cycle_cabinet | country / election_id / partyid),
                          data = dta_decades,
                          control=lmeControl(msVerbose = TRUE, sing.tol=1e-20))

# Save models because running the NLME models takes very long depending on the machine

save(reg_base_nlme, 
     reg_base_nlme2, 
     reg_pm_diss_nlme, 
     reg_single_party_gov_nlme, 
     model_decades_nlme, 
     file="regressions_nlme.RData")

### Create plots from NLME models ----


## Fig A2a: Effect plot: Basic model (NLME) ----
eff_basic_nlme <- data.frame(Effect("electoral_cycle_cabinet", reg_base_nlme,
                                    xlevels = 30))

ggplot(eff_basic_nlme, aes(x = electoral_cycle_cabinet,
                           y = fit, ymin = lower, ymax = upper)) + 
    geom_line() +
    geom_ribbon(fill = "grey50", alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change")


## Fig A2b: Effect plot: Single party government (NLME) -----
eff_single_party_gov_nlme <- data.frame(Effect(c("electoral_cycle_cabinet", "single_party_cabinet"), 
                                               reg_single_party_gov_nlme,
                                               xlevels = 30))

ggplot(eff_single_party_gov_nlme, aes(x = electoral_cycle_cabinet,
                                      y = fit, ymin = lower, ymax = upper, 
                                      fill = factor(single_party_cabinet), 
                                      linetype = factor(single_party_cabinet))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), 
                      name = "Single Party Government", 
                      labels = c("No", "Yes"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "Single Party Government", 
                          labels = c("No", "Yes"), 
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.4, 0.7))

## Fig A2c: Effect plot: PM dissolution power (NLME) ----
eff_pm_diss_nlme <- data.frame(Effect(c("electoral_cycle_cabinet", "diss_pm"), reg_pm_diss_nlme,
                                      xlevels=list(diss_pm=c(0,10), 
                                                   electoral_cycle_cabinet = seq(0,1,0.01))))

ggplot(eff_pm_diss_nlme, aes(x = electoral_cycle_cabinet,
                             y = fit, ymin = lower, ymax = upper, 
                             fill = factor(diss_pm), 
                             linetype = factor(diss_pm))) + 
    geom_line() +
    geom_ribbon(alpha = 0.3,
                aes(ymin = lower, ymax = upper)) + 
    scale_fill_manual(values = c("grey60", "dodgerblue2"), 
                      name = "PM Dissolution Power", 
                      labels = c("Lowest (0)", "Highest (10)"),
                      guide = guide_legend(reverse = TRUE)) +
    scale_linetype_manual(values = c(1, 2), 
                          name = "PM Dissolution Power", 
                          labels = c("Lowest (0)", "Highest (10)"),
                          guide = guide_legend(reverse = TRUE)) +
    labs(x = "Electoral Cycle", y = "Expected Absolute Poll Change") +
    theme(legend.justification = c(0, 0),
          legend.position = c(0.52, 0.7))


## Fig. A3: Historical development of the electoral cycle effect (NLME) ----

effect_decade_nlme <- as.data.frame(Effect(c("electoral_cycle_cabinet", "elecdecade"),
                                           model_decades_nlme,
                                           xlevels=list(elecdecade=unique(dta_decades$elecdecade), 
                                                        electoral_cycle_cabinet=seq(0,1,.01)))) %>% 
    mutate(elecdecade_print = paste("Decade:", elecdecade))

ggplot(data = effect_decade_nlme, aes(x = electoral_cycle_cabinet, y = fit)) +
    geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey80") +
    geom_line() + 
    facet_wrap(~elecdecade_print) +
    scale_x_continuous(breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1)) +
    xlab("Electoral Cycle") +
    ylab("Expected Absolute Poll Change") 

#### Create tables of NLME models ----


# get number of groups
groups_co <- length(unique(reg_base_nlme$groups$country))
groups_el <- length(unique(reg_base_nlme$groups$election_id))
groups_party <- length(unique(reg_base_nlme$groups$partyid))

## Table A5 ----

stargazer(reg_base_nlme, reg_base_nlme2, type = "text")

labels_nlme1 <- c("El. cycle", 
                  "El. cycle$^2$",
                  "El. cycle$^3$", 
                  "El. cycle", 
                  "El. cycle$^2$",
                  "GDP change", 
                  "Party support at last election",
                  "Election year - 1986",
                  "Single party gov.",
                  "El. cycle $\\times$ Single party gov.",
                  "El. cycle$^2$ $\\times$ Single party gov.",
                  "PM diss. power",
                  "El. cycle $\\times$ PM diss. power",
                  "El. cycle$^2$ $\\times$ PM diss. power",
                  "(Intercept)")

stagazer_nlme1 <- stargazer(reg_base_nlme, 
                            reg_base_nlme2, 
                            type = "latex",
                            dep.var.caption = "",
                            single.row = FALSE,
                            label = "reg:nlme1",
                            model.numbers = FALSE,
                            star.cutoffs = c(0.05, 0.01, 0.001),
                            column.labels = c("Model 5","Model 6"),
                            dep.var.labels.include = FALSE,
                            table.placement = "!h",
                            font.size = "scriptsize",
                            covariate.labels = labels_nlme1,
                            omit.stat = c("aic", "bic"),
                            title = paste0("Multilevel mixed-effects linear regression models of change in party support 
                                           (coefficients with standard errors in parentheses). The regressions replicate Models 1--3, but 
                                           take into consideration the autoregressive correlation structure (using the nlme R package).\\ 
                                           Groups (in all models): ", groups_co, " countries; ",
                                           groups_el, " elections; ", groups_party, " parties."))

#cat(stagazer_nlme1, sep = "\n", file = "regression_nlme1.tex")

## Table A6 ----

stargazer(reg_single_party_gov_nlme, 
          reg_pm_diss_nlme, type = "text")


labels_nlme2 <- c(
    "El. cycle", 
    "El. cycle$^2$",
    "GDP change", 
    "Party support at last election",
    "Election year - 1986",
    "Single party gov.",
    "El. cycle $\\times$ Single party gov.",
    "El. cycle$^2$ $\\times$ Single party gov.",
    "PM diss. power",
    "El. cycle $\\times$ PM diss. power",
    "El. cycle$^2$ $\\times$ PM diss. power",
    "(Intercept)")

stagazer_nlme2 <- stargazer(reg_single_party_gov_nlme, 
                            reg_pm_diss_nlme,
                            type = "latex",
                            dep.var.caption = "",
                            single.row = FALSE,
                            label = "reg:nlme2",
                            model.numbers = FALSE,
                            star.cutoffs = c(0.05, 0.01, 0.001),
                            column.labels = c("Model 7", "Model 8"),
                            dep.var.labels.include = FALSE,
                            table.placement = "!h",
                            font.size = "scriptsize",
                            covariate.labels = labels_nlme2,
                            omit.stat = c("aic", "bic"),
                            title = paste0("Multilevel mixed-effects linear regression models of change in party support 
                                           (coefficients with standard errors in parentheses). The regressions replicate Models 1--3, but 
                                           take into consideration the autoregressive correlation structure (using the nlme R package).\\ 
                                           Groups (in all models): ", groups_co, " countries; ",
                                           groups_el, " elections; ", groups_party, " parties."))

#cat(stagazer_nlme2, sep = "\n", file = "regression_nlme2.tex")

## Table A7 ----
labels_decades_nlme <- c("El. cycle", "El. cycle$^2$",
                         "El cycle$^3$",  
                         "Decade – 1970s", "Decade – 1980s",
                         "Decade – 1990s", "Decade – 2000s", 
                         "Decade – 2010s",
                         "GDP change",
                         "Party support at last election",
                         "El. cycle $\\times$ Decade – 1970s.",
                         "El. cycle$^2$ $\\times$ Decade – 1970s.",
                         "El. cycle$^3$ $\\times$ Decade – 1970s.",
                         "El. cycle $\\times$ Decade – 1980s.",
                         "El. cycle$^2$ $\\times$ Decade – 1980s.",
                         "El. cycle$^3$ $\\times$ Decade – 1980s.",
                         "El. cycle $\\times$ Decade – 1990s.",
                         "El. cycle$^2$ $\\times$ Decade – 1990s.",
                         "El. cycle$^3$ $\\times$ Decade – 1990s.",
                         "El. cycle $\\times$ Decade – 2000s",
                         "El. cycle$^2$ $\\times$ Decade – 2000s",
                         "El. cycle$^3$ $\\times$ Decade – 2000s",
                         "El. cycle $\\times$ Decade – 2010s",
                         "El. cycle$^2$ $\\times$ Decade – 2010s",
                         "El. cycle$^3$ $\\times$ Decade – 2010s",
                         "(Intercept)")


# get number of groups
groups_co_dec <- length(unique(model_decades_nlme$groups$country))
groups_el_dec <- length(unique(model_decades_nlme$groups$election_id))
groups_party_dec <- length(unique(model_decades_nlme$groups$partyid))

stargazer(model_decades_nlme, type = "text")

stagazer_nlme3 <- stargazer(model_decades_nlme,
                            type = "latex",
                            dep.var.caption = "",
                            single.row = TRUE,
                            label = "reg:nlme3",
                            model.numbers = FALSE,
                            star.cutoffs = c(0.05, 0.01, 0.001),
                            column.labels = c("Model 9"),
                            dep.var.labels.include = FALSE,
                            table.placement = "!h",
                            font.size = "scriptsize",
                            covariate.labels = labels_decades_nlme,
                            omit.stat = c("aic", "bic"),
                            title = paste("Multilevel mixed-effects linear regression models of change in party support 
                                          (coefficients with standard errors in parentheses). The regressions replicates Models 4, but 
                                          takes into consideration the autoregressive correlation structure (using the nlme R package).\\ 
                                          Groups: ", groups_co_dec, " countries; ",
                                          groups_el_dec, " elections; ", groups_party_dec, " parties."
                                          , sep = ""))

#cat(stagazer_nlme3, sep = "\n", file = "regression_nlme3.tex")
